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We develop an approximate formalism suitable for performing simulations of the thermal dynamics 
of interacting Bose gases. The method is based on the observation that when the lowest energy 
modes of the Bose field operator are highly occupied, they may be treated classically to a good 
approximation. We derive a finite temperature Gross-Pitaevskii equation for these modes which is 
coupled to an effective reservoir described by quantum kinetic theory. We discuss each of the terms 
that arise in this Gross-Pitaevskii equation, and their relevance to experimental systems. We then 
describe a simpler projected Gross-Pitaevskii equation that may be useful in simulating thermal 
Bose condensates. This classical method could be applied to other Bose fields. 
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I. INTRODUCTION 



CN , The achievement of Bose-Einstein condensation (BEC) in a dilute gas [Q-g[ offers the possibility of studying the 

^ ' dynamics of a quantum field in the laboratory. In principle it presents an opportunity to directly compare computa- 
^^ tional predictions with experimental results; however carrying out dynamical calculations of thermal quantum fields 
, -- is an extremely difficult problem that generally requires severe approximations to be made. 

r The most successful finite temperature theories of BEC are based on second-order perturbation theory, and are 

f~^ limited to the case of thermal equilibrium away from the region of critical fluctuations [y-|6| . These theories have 
T-H allowed the accurate determination from first principles of quantities such as excitation frequencies and damping rates 

of Bose-condensed systems. However, their extension to dynamical situations is computationally difficult. 

At very low temperatures when most of the atoms are in the condensate, the Gross-Pitaevskii equation (GPE) has 

proved remarkably successful in numerically modelling BEC experiments. The time-dependent GPE has the form 



^ ■ '^dT " ^^P'/'W + f/o|^(x)PV;(x), (1) 



^ ' where f/g = 47r?i a/m is the effective interaction strength at low momenta, a is the s-wave scattering length, and m 
is the particle mass. Hsp is the single particle Hamiltonian 



h^ 



^. Hs, = - — V' + V,,,^i^), (2) 



where Vtrap(x) is the external confining potential of the system. The GPE can be derived by a number of different 
approaches (e.g. a number-conserving approach 0), but a direct method is to take the mean mean value of the 
equation of motion for the Bose field operator and assume that the quantum fluctuations can be neglected. This 
procedure effectively assumes that the field is well approximated by a coherent state. 

It has been argued that the GPE can also describe the dynamics of a Bose-Einstein condensate at finite temperature 
[p[-pd| for the reason that in the limit that all the low- lying modes of the system are highly occupied {Nk ^1), the 
classical fluctuations of the field operator $(x, t) will be much larger than the quantum fluctuations. It is therefore 
reasonable to neglect the quantum fluctuations, and thus all highly-occupied modes may be well approximated by a 
coherent wave function. This is analogous to the situation in laser physics where the highly occupied laser modes 
can be well described by classical equations. Despite this proposal appearing in the literature in 1991 [p[, very few 
numerical calculations have been performed. The first was by Damle et al. [ |l2| , who calculated the approach to 
equilibrium of a near-ideal superfluid using the GPE. Subsequently, Marshall et al. jl^ performed two-dimensional 
simulations of evaporative cooling of a thermal Bose field in a trap. More recently, papers by Stoof and Bijlsma |lj] 
and Sinatra et al. Pq | have used classical methods in dynamical calculations of thermal one-dimensional Bose-Einstein 
condensates. Coral et al. |1^ have performed dynamical calculations that treat several modes of a 3D homogeneous 
Bose gas classically. Their method, while not specifically employing the GPE, is similar to the approach we suggest 
here and for which we have presented our own numerical results in reference p^ . Classical approximations to other 

1 



quantum field equations have also been successful in the calculation of the dynamics of the electroweak phase transition 

@- ^ 

The major advantage of using the GPE to describe thermal dynamics is simply that, while it is still a major 
computational task, it is possible to numerically solve the equation for realistic systems in a reasonable amount of 
time. In addition, the GPE is non-pcrturbative and it should be possible to study the region of the BEG phase 
transition, where perturbation theory often fails. 

There are, however, a number of problems associated with the use of the GPE to represent the entire Bose field at 
finite temperature. It is a classical equation, and so in equilibrium it will satisfy the equipartition theorem — all modes 
of the system will have an occupation of N^ = ksT/ck- Thus, if we couple the GPE to a heat bath and numerically 
solve the equation with infinite accuracy, we will observe an ultra-violet catastrophe. Also, we can see that the higher 
the energy of any given mode, the lower its occupation will be in equilibrium and at a sufficiently high energy the 
criterion Nk ^ 1 will no longer be satisfied. For these low occupation modes a form of kinetic equation is more 
appropriate. The solution to both of these problems is to introduce a cutojf in the modes represented by the GPE. 

In this paper we develop an approximate formalism in which the low-lying modes of the system are described 
non-perturbatively by the GPE, coupled to a thermal bath described by a quantum Boltzmann equation. We derive 
a finite temperature Gross-Pitaevskii equation (FTGPE), and discuss the terms that couple the part of the field 
operator represented by a coherent wave function to the thermal bath. In particular we show how a description of 
loss via elastic collisions arises naturally in the formalism. 

Several other authors have developed formalisms for non-equilibrium dynamics using quite different theoretical 
methods — we mention Gardiner and ZoUer [p^-pl[, Proukakis et al. |23], Stoof |23 24|, Zaremba et al. pa], Walser 
et al. pq] , and finally Sinatra et al. p7[ . The work we present here has elements in common with several of these. 
In particular, the average of the quantum Lan gev in equation written down in the formalism of Stoof p3,p3] would 
correspond to the FTGPE we derive in section IV . 

This paper is organized as follows. In sectionO we write down and discuss the Hamiltonian which is our starting 
point. In section HI we outline how we decompose the field operator into a coherent and incoherent region, and 
in section |V| we derive a finite temperature Gross-Pitaevskii equation that forms the main result of this paper. In 
section M wc discuss the terms that arise in the FTGPE, their relation to experiments, and their approximate form 
in terms of occupations numbers of incoherent region modes. We discu ss a simple finite temperature equation which 
we call the projected GPE in section VI and finally conclude in section VIE 



II. HAMILTONIAN 



We begin with the usual second quantized many-body Hamiltonian for a system of identical, structureless bosons 
with pair-wise interactions 



H — Ho + Hj, 



Hi = ^ f d^x /rf3x',|,t(x^i)^|ft(x',t)F(x-x')^(x',i)^(x,i). 



(3) 

(4) 

(5) 



The non- interacting part of the Hamiltonian, Hq, corresponds to an ideal gas and can be diagonalized exactly by the 
eigenvectors of Hsp. The quantity Hi describes two-body interactions via the interatomic potential V^(x). The field 
operator ^(x, t) annihilates a single boson of mass m at position x and time t, and obeys the equal-time commutation 
relations 



*(x,i),^(x',t)j = [*^(x,i),^t(x',t)l =0 



(6) 



*(x,i),'l't(x',i) =,5(x-x' 
The field operator is normalised such that 

/a*.(.,OT(x,o^a', 

where N is the particle number operator of the system. 



(7) 



(8) 



A. Basis set representation 

It is useful to expand the field operator on a basis set 



(9) 



where (^„(x) is a mode function, and a„(i) annihilates a particle in mode n at time t. These operators obey the 
equal-time commutation relations 



[am, an] = [ajjjjay = 0, 



(10) 



[^m? ^nj ^^ ^f 



H = ^hujnaldn + ^ X! (P<l\V\mn)alalarnan, 



pqmn 



where we have defined the symmetrized matrix element 



1 



{pq\V\mn) - - / d^x / d^x' 0;(x)0;(x')F(x-x')0™(x')0„(x) 
+ i /d^x /dV,/)!(x)0*(x')F(x-x')0„(x')</>™(x). 



(11) 



where we have dropped the time labels for clarity. If we substitute equation (y) into the Hamiltonian 1^) and take 
the set {(/)„} to be the eigenvectors of Hgp, we find 



(12) 



(13) 



whose use significantly reduces the length of the equations of motion we later derive. Equation (13) represents both 
direct and exchange collisions, which are physically indistinguishable for identical bosons. 
The Heisenberg equation of motion for the individual mode operator Op is therefore 



da 



h—^ = hajpUp + y^ {pq\V\mn)alaman, 



(14) 



and we can define slowly-varying operators 



Cir) (Ji'nfi 



Auj^t 



SO that the equation of motion for the annihilation operator becomes 



da 



i^—rf- = y^^{pq\V\mn)alaman( 



,i(uJp-\-UJq—UJj-ri—'-^7z)i 



q7nn 



(15) 



(16) 



B. Effective low-energy Hamiltonian 



The Hamiltonian described above contains spatial integrals over the bare interatomic potential between two atoms, 
V^(x). However, it is well-known that at low temperatures the scattering of neutral atoms in three dimensions can be 
described by the s-wave scattering length a. This parameter is often introduced into the theory by replacing the real 
interatomic potential by the contact potential 



l/(x-x')^;7o(5(x-x'), Uo 



Anh a 



(17) 



The interaction strength Uq can be shown to arise from the increase in kinetic energy of a two-particle wave function, 
when an excluded region of radius a is introduced corresponding to a hard sphere interaction potential pq |. The 
contact potential approximation, however, can lead to ultraviolet divergences in theories of BEC if it is simply 



substituted into the Hamiltonian (H). This is not surprising, as the delta-function potential can scatter high-energy 
atoms just as effectively as low-energy atoms. Physically this is unrealistic, as momentum transfer between atoms 
will vanish at high momenta (k > 1/a). The contact potential is a low-energy approximation, and care must be taken 
when summing over high energy states. 

The ultraviolet renormalization of the theory can be achieved by introducing the two-body T-matrix into the 
Hamiltonian, resulting in a high- momentum cutoff K to the the states considered yj. This procedure is valid as long 
as the cond ition Ka <C 1 is satisfied for the entire gas. We will consider the issue of ultraviolet divergence further in 



section VA 



III. PROJECTION OPERATOR 

The aim of this paper is to represent the highly occupied modes of the field operator ^(x) by a wave function ^(x). 
It is reasonable to neglect the quantum fluctuations of these modes, and therefore ■(/'(x) represents a classical field. 
To achieve our goal we divide our representation of the field operator into two separate regions. The first, in which 
the condition Nk ^ 1 is satisfied, we call the coherent region C. The other which contains the remainder of the field, 
is denoted the incoherent region /. 

We define the coherent region projection operator 



V=Y.W)W\, (18) 



vec 

where the region C is determined by the requirement that {a\,aiy) 3> 1, and the set \v) defines some basis in which 
the field operator is approximately diagonal near the cutoff. This condition is imposed simply so that in equilibrium 
the quantity {a]^ay) has a well-defined average at the boundary of C that can be interpreted as a mode occupation 
number. The position of the cutoff is a choice that must be made before any calculation, and is required for the 
construction of the initial wave function.^ 
Operating on the field operator with V gives 



7^*(x) = Y. '^-W Z^''^' C(x')^(x'), 



eeV3(x). (19) 

such that '0(x) is the field operator for the coherent region. We now introduce the orthogonal projector Q ~ I — V 
and define 

Q*(x) = ^afc(/)fe(x), 

s 7?(x). (20) 

The quantity 77(x) is the field operator for the incoherent region and represents an effective thermal bath. Quantum 
fiuctuations are important for these modes — in fact we will later assume that (0^) « for the large majority of k ^ C. 
The full field operator is 

^(x) = [•P + Q]*(x), 

= V'(x)+77(x), 

= ^ai.0^(x) + ^afc(/ifc(x), (21) 

where we indicate indices within C by Greek subscripts, and outside C by Roman subscripts. We shall follow this 
convention throughout the body of this paper. 

As an example of the size of the regions we consider a ^^Rb gas in a harmonic trap with a geometric mean trap 
frequency oi uj = 2tt x 100 Hz. For a condensate of 10^ atoms at a temperature of 640 nK (and hence a total number 
of trapped atoms of about 5 x 10^), we find that a quantum level with energy e — /i « 15huj has a mean occupation 
of (N) « 10, so this would be an appropriate boundary for the coherent region. In comparison, the remainder of the 
gas spans energies up to about E = IbOOtiuj, meaning that the incoherent region contains many more quantum states 
and is much larger than the coherent region. 



IV. EQUATIONS OF MOTION 

A. Hamiltonian 

We now substitute the decomposition of equation (|T]) into the Hamiltonian (||) . We assume ior k ^ C that (j)k (x) 
is an eigenstate of -ffsp, and so Hq simphfies to 

Hq = ^{a\Hsp\l3)al,ai3 + 'h^ujkalak. (22) 

a/3 k 

For the interaction part of the Hamiltonian we have 

Hi ^ -Y^ {a(3\V\x(T)alala^a^ (23a) 

a/3xcr 

+ Yl \^{af3\V\xn)alala^an + h.c (23b) 

+ o X! {(^P\V\'fnn)"'i^l"'man + h.c (23c) 



af3mn 



+ 2 Y {aj\V\xn)aid^a]a„ (23d) 

ajxn 

+ 2^ \{aj\V\mn)al^a'^^aman + h.c. (23e) 

ajmn 

+ - ^ (A:j|y|TOn)aj^a]a„o„, (23f) 

kjmn 

where the symmetrized matrix element {kj\V\mn) is defined in equation dl3), and h.c. stands for hermitian conjugate. 
Using ( |22|) and ( |23| ) we now derive the Heisenberg equations of motion for the operators in each region. 

B. Coherent region 

The evolution of a mode of the coherent part of the field operator is given by 

*?J^=^(a|i?.p|/?)a/3 (24a) 



dt 





Y{al3\V\x<j)dld^d„ (24b) 



Px'y 



+ ^(ag|y|xCT)4axa^ (24c) 

+ 2 ^ {a(3\V\ma)aldmaa (24d) 

(3ma 

+ y ^{al3\V\m.n)aQar^^an (24e) 

I3mn 

+ 2 ^ {aq\V\ma)d\d^d„ (24f) 

+ ^(ag|y|mn)aja,„a„. (24g) 

qnin 

in which the coupling to the incoherent region is made explicit. We now begin to introduce our approximations. The 
condition that a mode is in the coherent region C is that the population (d'ld^) ^ 1, which allows us to neglect the 
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FIG. 1. Representation of the degree of coherence in the different regions of the Bose field. The shading indicates schemat- 
ically the coherence of the field. The classical region of the field in indicated by C, and the incoherent region by 7. The states 
in / near the boundary of C will be partially coherent. 



quantum fluctuations of the projected field operator V'(x). Thus we assume that the region C is well approximated 
by a mode function given by the mean value 



V'(x) = (^(x)), 
and we can expand the wave function on our basis functions as 

^(x) = ^(a^)(/i^(x) 

= ^Ci.0,.(x). 



(25) 



(26) 



vec 



We obtain the finite temperature GPE by taking the mean value of equation (|2J) . In this procedure we expand the 
coherent region operators as a,y = c,^ + S^, and then neglect the terms involving the quantum fluctuations S^. A 
typical term simplifies as follows 



9XCT 



9XO- 



X^cr 



(27) 



where the matrix element is time dependent as the wave function ip is not a stationary state. 

The incoherent region / is, for the most part, best represented by number states. However, this is not always 
the case. In particular the states within / but near the boundary of the two regions will be partially coherent, as is 
illustrated in figure |l|. The expectation value (aj) in this transitional region will not be zero, and so terms such as (|27| ) 
are retained in our equations. This is different from other mean field theories in which all coherence is absorbed into 
the GPE. In systems that are partially condensed, however, the effect of these terms will be small, as the transition 
region will be narrow compared to the full width of the incoherent region. 

The full equation of motion for a coherent region mode obtained by taking the mean value of all terms in equa- 
tion iM) is 



dCa 






(28a) 
(28b) 
(28c) 



+ 2^(ai/;|y|mi^)(a„) 



(28d) 



+ ^{ai^\V\mn){a,nan) (28e) 

+ 2^(ag|l/|mV')(4a™) (28f) 

qm 

+ '^{aq\V\mn){a\aman)- (28g) 

qmn 

We can convert this expression to the spatial representation by applying the operation VLg(^ |a) to both sides. Using 
the contact potential approximation and recognizing X^qpc I'-'^)('^I ^^ '^^^ projector of (|l8[), this procedure results in 
an equation we call the finite temperature Gross-Pitaevskii equation (FTGPE) 

i^^^ = HspH^) + UoP {|^(x)pV;(x)} (29a) 

+ UoV {2|V(x)P(??(x)) + V'(x)2(,)t(x))} (29b) 

+ UoP {7/^*(x)(7}(x)7}(x)} + 2V^(x)(77t(x)7)(x))} (29c) 

+ C/oP{(f/t(x)7}(x)77(x))}, (29d) 

where f}(x) is defined by (EQ). The FTGPE constitutes the main result of this work, and the remainder of this paper 
is devoted to discussing the physics of this equation. 

The only approximation that has been made in the derivation of the FTGPE is that the modes it represents must 
be highly occupied. No perturbative techniques have been used, and therefore the equation should be valid as long 
as the condition Nk ^ 1 is satisfied. As a corollary, we expect the FTGPE could be used to study the region of the 
phase transition in certain circumstances. 

The FTGPE describes the full dynamics of the region C and its coupling to an effective thermal bath 77(x). The 
initial wave function for the coherent region will be made up of a sum over a basis with amplitudes with random 
phases as is appropriate for a thermal system lO]. Despite the fact that the FTGPE is completely unitary and 
reversible, we expect that it will evolve general initial states of V'(x) to an equilibrium determined by the temperature 
and chemical potential of the field 7?(x). This is because deterministic nonlinear systems exhibit chaotic, and hence 
ergodic, behavior if more than a few degrees of freedom are present p9| . In fact, we have shown in reference p7| that 
a simplified form of the FTGPE that we discuss in section W^ can indeed describe evolution towards equilibrium. 

Each of the lines of the FTGPE represents collision processes involving a different number of coherent region states. 
We describe them briefly here and then discuss them in more detail in section |V|. 

The terms on the first line of the FTGPE (p9q) represent purely coherent region dynamics. The first term describes 
the free evolution of the wave function ip, while the second represents evolution due to two particles from C colliding, 
with both particles remaining inside the coherent regi on. 



The terms on the second line of the FTGPE (29b) we refer to as the linear terms. These describe two coherent 
atoms interacting, resulting in one remaining in C and one escaping to the incoherent region (and the reverse process) 
and are depicted in figure 0(a). These are stimulated processes as the terms contain three coherent region labels, and 
they result in the transfer of some coherence to the bath ?7(x) (see figure |l|). However, because the coherent region is 
much smaller than the incoherent region these terms can often be neglected in comparison with the third and fourth 
lines. 



The third line of the FTGPE (29c) will usually be more important than the second. The first term, which we 
call the anomalous term, represents the collision of two coherent atoms with two incoherent atoms resulting as in 
figure M(b) . If the region C represents only a single condensate in thermal equilibrium then this term cannot conserve 
energy and therefore it cannot describe real processes]^ However, it can describe uirtua/ processes and thus contributes 
to the appearance of both the two-body and many-bo dy T -matrices. 

The second term of the third line of the FTGPE ( ^9c| ) represents a coherent atom colliding with an incoherent 
atom, with one atom remaining in each region after the interaction as in figure 0(c). In equilibrium this process 



represents the mean field of the incoherent region acting on C, and can be added to the |^(x)p term of (29a). Away 
from equilibrium this term also describes scattering processes, identical to those described in the model of condensate 
growth developed in references P^'HI ■ 



'^When the coherent reg ion is made up of two or more condensates then the first term of (29c) can describe real processes as 



we discuss in section VB 
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FIG. 2. An illustrative diagram of the collision processes represented by the terms that appear in the finite temperature 
GPE ( pgj ) that contain r\ and/or r]^ . (a) The terms linear in r\. (b) The anomalous term, (c) The scattering term, (d) The 
growth term. The labels p and q indicate initial states, and the labels m and n are scattered states. The reverse processes are 
also possible. 



Finally the fourth line (29d) represents the collision of two incoherent atoms in which one is transferred to the 
coherent region C as depicted in figure g(d) . This is the growth process described in references ||3yj3l[ , and is the 
main contribution to the transfer of population between the coherent and incoherent regions. 



C. Incoherent region 

The coherence of the majority of levels outside of C is negligible, and therefore most of the incoherent region is 
approximately diagonal in the number state representation. The energy of the quantum levels is large enough that 
the mean field of the wave function ?/; does not significantly affect this region, and so we assume that Hi is a small 
perturbation to Hq. Therefore quantum kinetic theory can accurately describe the evolution of the majority of this 
part of the quantum field, with the appropriate modifications to treat the coupling to the coherent region. 

We can derive an equation of motion for the incoherent region using similar techniques to those used in the 
derivation of the quantum Boltzmann equation (QBE) and employing two major approximations — the random phase 
approximation (RPA) and Wick's theorem [Q . We demonstrate this approach in appendix ^ where we derive QBE 
from the field theory equations. We also use these methods to simplify the thermal terms of the FTGPE as we discuss 
in the next section. The kinetic equation for the incoherent region is not the main topic of this paper; however, for 
completeness we write it down and discuss the terms that arise in appendix O. 

V. DETAILED ANALYSIS OF THE COUPLING TERMS OF THE FTGPE 



In this section we interpret in detail the meaning of, and find expressions for, the terms involving mean values of 
combinations of the bath operator 77(x) on the right-hand side of the FTGPE, equation (p9|). As the incoherent region 
is best represented by the number occupation of the quantum levels, it is appropriate to express these mean values in 
terms of the occupation numbers and to do so we are guided by kinetic theory. 



Expressions are required for the quantities (f}), (fj^), {fiTJ), {v fj), and {fyfjfj). These can be derived from the equation 
of motion for the corresponding combinations of creation and annihilation operators for the incoherent region. The 
equations we begin with are rather long, and we write them out fully in appendix y. 

We then follow the general procedure used to derive the QBE (see appendix |A|) . First we approximately integrate 
the starting equations for the combinations of creation and annihilation operators, before taking the mean value 
and employing Wick's theorem and the RPA to simplify the resulting terms such that they depend explicitly on 
populations Up — {alap). These expressions are then substituted back into (^), {if), etc, and then into the FTGPE. 
Each of the terms requires slightly different treatment, and we provide key details in the following subsections. 

A. The linear terms 

The two terms involving {fj{x)) or (77''^(x)) in equation (poj ) describe the collision of two coherent atoms, where 
one particle remains in the region C and the other is transferred to the incoherent region (along with the reverse 
process). In systems where there is significant population in the incoherent region, these terms will not be very large 
in comparison to the terms of higher order in r) due to the requirements of energy conservation and the relatively 
small size of the region C. 

Beginning with equation ( pl| ), eliminating the free evolution via the transformation fip — ape'''^p*, taking the mean 
value, and using Wick's theorem and the RPA, the only terms that survive are 



d{ap) 
dt 



= (pVI^IV'V')e"^''' + 2^(pq|F|gV)n?e*'^''*- (30) 



We can also drop the last term in this equation using the rotating wave approximation, as the components of -0 will 
be oscillating at frequencies smaller than LOp. Thus the result is 

ifi ^^ = {p^\V\ii^^P)e'''^K (31) 

In order to find a approximate solution for equation ( pl| ) , we assume that we can expand the coherent region wave 
function in a basis that is approximately diagonal — essentially a quasiparticle basis. We write 

V'(x) = ^g.C.(x)e— "*, (32) 

(7 

where the coefficients c^ ~ Co-e*'^'^*. The equation of motion for dp becomes 

*^^ = E^^2^^.(p/3|y|xa)e--^-*, (33) 

where we have introduced the notation 

LOpqmn = ^p + ^q ^ ^m ~ tOn- (34) 

If the set I^ctI is a good basis for the coherent region, then the exponential term contains most of the time dependence 
of equation (|33|). We can therefore take everything else outside the integral, and use the standard result 

^ p(-] -iirSix), (35) 



X + le \x 

to find an approximate solution to equation (p3|). Incorporating the free evolution in the solution, we find 

(op) = 2^ < ^ iT^Cf-jC^CcripfilVlxajSihujpfj^a) > ■ (36) 

/3xo" PPX J 

Since we are mainly interested in the kinetic processes that can occur, we neglect the energy shift de scrib ed by the 



principal part (the first term on the right-hand side). The relevant term in the FTGPE in basis form (28d) is thus 



dc 
«^-^ = ••• -27ri^(aK|T/|pi/)c*c^^4cxC^(p/3|V|xcr)^(?iWp/3;^^), (37) 

where we have expanded all the condensate wave functions as in (JS^)- We note that equation (|3^) contains only 
coherent region amplitudes {ca} because the processes we have included are all stimulated. 

A situation where the terms discussed in this section may become important is for experiments in which a Bose 
condensate near T = is disturbed by a sudden change in its scattering length by the use of a Feshbach resonance. 
Such experiments have been carried out by Donley et al. using ®^Rb at JILA ||3^. The change in scattering length 
causes the collision processes represented by the linear terms of the FTGPE to become energetically allowed. This 
description of the physics is identical to the argument by Duine and Stoof that the loss observed can be explained 
via elastic collisions due to an imaginary part of the many-body T- matrix |Q . We will address this issue further in 
a future paper. 

B. The anomalous term 



We now consider the term involving (777}) on the third line of the FTGPE (29c). Expanding this quantity in the 
incoherent region basis gives 

(W) = Y. 0m0„(ama„)e-'(""+'^")*. (38) 

mn 

To find an expression for (dmO'n) we use the equation of motion for this quantity given as equation ( |C2| ). Eliminating 
the free evolution, taking the mean value, expanding the coherent region wave function in the approximately diagonal 
basis according to equation (|3^), and finally using Wick's theorem and the RPA we obtain 

+ ^(mn|y|fcj)(&feaj)e'"™"''^*(l+n,„ + n„). (39) 

kj 

The first line of this equation describes two particles scattering from the coherent region into states m and n. The 
second line would usually be ignored in the RPA, as it is of higher order than the first line. However the matrix element 
of this line describes particles from (rn,n) scattering into {k,j) and then onto other states. This offers the possibility 
of ladder diagrams, such that the two particles could scatter back into the coherent region without interacting with a 
third atom. We retain this term in equation (p3) as the mean value it contains is of the same form as the that on the 
left-hand side. Because of this, equation ( |39| ) has the form of a Lippmann-Schwinger equation in the time domain, 
and we will see that this is where the T-matrix appears in the formalism. 

To solve this equation we start by considering the first line only, as it is the lowest order term. Once again, most 
of the time dependence is contained within the exponential, and so we find the solution is 

,. , , _ -e;^ . , {mn\V\x(T)il + n„ + n„)e'"'"">^'"* 



^CT ^X ^ ^o- ^m c„ 



-i7r^c^Ccr(mn|t/|xo-)(5(ex + e„ - tm - Cn)- (40) 

For a single condensate near thermal equilibrium, the energy delta-function can never be satisfied as it requires two 
low-energy, coherent atoms from within the coherent region to collide and result in two high-energy, incoherent atoms. 
We will return to this point later in the section. We assume that the full solution of equation ( p9| ) is of the same form 
as (B^) but with the interaction potential V replaced by a T-matrix T"p 

(amfln) = 2^ C^Ca . (41) 



This is a solution of equation (^9|) if the operator T^p obeys 

1 + Uk + Uj 



T^^{z)^V + Y,V\kj) 3 ''[_ ■' (fcj|r"P(z), (42) 

kj 
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where we have identified the parameter z = e^^ + ta + iS ds the incoming energy of the two particles in the coUision. The 
small imaginary part id in this parameter generates the delta function term in ([40|). Equation (|42| ) is the definition 
of a restricted many-body T-matrix as the indices k,i are defined to be outside the coherent region. 

The many-body T-matrix describes collisions in the presence of a medium. It takes into account the fact that the 
virtual states that two particles pass through in the collision may be occupied, and the scattering rate enhanced by 
a factor of (1 + rife + rij). However, as the states fc, j are in the incoherent region the populations are generally small, 
and if we can approximate Uk = rij = we recover the definition of a restricted two-body T-matrix. 

If we now substitute the solution equation ([4l|) into the basis set FTGPE equation (Eq), we find from the terms 



(28t) and (280 



^,dc^ _ , Y-„ „ /„./.,T.u,_^ . Y-/„.,.,T.,™„vY-„ „ {mn\T-P\xa) 



i^-if = ■ ■ ■ + ^CxC^{a^p\V\x(j) +^{aip\V\mn)^CxCa 



dt 

(77in|T"P 



Xcr mn xa ^ 



^CxC<,(aV'| 



V ^^V\mn)- 



\xo). 



= (aV'|T"P|V'V')- (43) 

Thus the anomalous term introduces the restricted T-matrix into coherent region collisions, and it is appropriate to 
approximate this T-matrix by a contact potential y . With this treatment of the anomalous term we have carried 
out the ultra-violet renormalization of the theory. 

It is useful to make two additional points. First, the T-matrix that enters our equations is not the full two-body 
T-matrix as the sum over intermediate states only includes levels in the incoherent region. Thus the contact potential 
we use in our calculations should be 

r"P^C7o^(x-x'), (44) 

where Uq ^Vq. The remainder of the terms that would 'upgrade' C/q to Uq are included directly in the simulation of 
the coherent region using the FTGPE, and so all T-matrix effects are actually included in the formalism. 
In practice wc find these two quantites are related in the homogeneous limit H] by 



where aK is defined by 



^ - raK 



and K is the wave vector of the cutoff between the coherent and incoherent region. Thus as long as the condition 
UoaK ^ 1, or equivalently Ka <C 1 is satisfied, we are justified in setting C/q ~ C/o. As in any calculation the coherent 
region will be rather small, and experimentally measured scattering lengths for alkali atoms are generally not known 
with an accuracy of better than 10%, this approximation is well justified in the homogeneous case. It seems reasonable 
to expect the same condition should hold true in a trap. 

The second point is that the remaining terms of the FTGPE still have the interatomic potential rather than the 
T-matrix in their matrix elements. It turns out that it is reasonable to use the contact potential approximation in 
these, although we have not explicitly upgraded the matrix elements to T-matrices. For further details we refer the 
reader to references ||35|,p|, p2| , p6[ 

1. Elastic loss in condensate collisions 



In section |VB| we stated that the delta function in the solution of the anomalous term (^0|) could not be satisfied for 
a system with only one condensate near thermal equilibrium. The situation is different when we consider the collision 
of two condensates, where the T-matrix can have an imaginary part. 

In the formalism described here the coherent region C is defined such that it contains only the modes of the gas 
whose occupation number satisfies N^ ^1. If we consider a condensate that has been formed in a harmonic trap, 
but then quickly released into free space, we can analyse the wave function in terms of a plane-wave basis. The region 
C will typically be defined by a small spherical or ellipsoidal region in fc-space about a central wave-vector K. For 
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FIG. 3. A depiction in fc-space of two condensates colliding in the centre-of-mass frame. The two shaded areas centred about 
±K/2 indicate the coherent region C. The arrows indicate a possible collision process between two coherent particles, in which 
both end up outside the coherent region. The dashed circle indicates the region of all possible collisions that conserve both 
energy and momentum. 

a condensate released at rest we have K = 0, but by applying a Bragg-pulse to the condensate on release, it can 
be split into two equal parts — one at rest and one with a momentum KK. ||37| . In the centre-of-mass frame, the two 
condensates have momenta —KK./2 and +TiK./2 respectively. 

Now consider this system analysed in the plane-wave basis. If KK. is large compared to the momentum width of 
the two condensates, then the region C will be made up of two distinct parts of fc-space, as depicted in figure |[ This 
means that in the collision of the two condensates it is possible for an atom from each to collide, and then scatter 
into any spatial direction while energy is still conserved. A large number of these collisions will take both particles 
outside the region C, as depicted in figure ^ As these processes cause scattering of bothparticles into modes that are 
otherwise empty (and are hence spontaneous), they cannot be represented by the GPE (|I|) (see appendix^ for further 
explanation). However, these collisions can be represented in the FTGPE via the anomalous term {fji)). Unlike the 
case of a single condensate, the delta function part of equation (^2) can now be satisfied, and therefore real transitions 
can occur. 

We can also see that if the relative momenta of the two condensates is not large, then the two parts of the region 
C will overlap and most of the circumference of the dashed circle in figure p^ will lie within C . This means that the 
GPE can describe condensate collisions at low momenta for which spontaneous collisions can be neglected. However, 
at high relative momenta other methods, such as described here, are required. 

The process of elastic loss is a form of spontaneous Beliaev damping. An analogous phenomenon for a trapped 
BEG is, for example, when a high energy coherent collective excitation is generated in a ground state condensate. 
The coherent region is again divided into two parts, and the excitation can interact with the ground state condensate 
and spontaneously decay into two lower energy quasiparticles. 

2. Four-wave mixing 

In an elegant experiment, the Phillips group at NIST demonstrated the atom-optical equivalent of four-wave mixing 
with Bose-Einstein condensates ^^l- After releasing a condensate from a trap, two separate Bragg pulses were applied 
in succession such that the BEG split into three distinct parts, each with a different momentum (two moving and one 
at rest in the lab frame). For appropriately chosen momentum values, a fourth condensate was observed to appear. 

This experiment can be understood by considering figure pi but now with a third condensate at the tip of the arrow 
at the top of the dashed circle. While the entire dashed region is still available in collisions between atoms from the 
first two condensates, the presence of third condensate stimulates transitions into this particular mode, resulting in 
the formation of a condensate at the tip of the downward pointing arrow. 

The ordinary GPE (|l|) can describe the formation of the fourth condensate as this is a stimulated collision process 
|p9^,Eo[. However, the new condensate that appears is entangled with the atoms that are scattered into the third 
condensate as the colliding atoms are necessarily correlated. This correlation cannot be described by the ordinary 
GPE, and other methods are required pi. We emphasize that such effects can be described by the anomalous term 
in the FTGPE. 
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C. The scattering term 

We now consider the term containing {rff}) in the third line of the FTGPE (p9q). To find an explicit expression for 
this quantity, we first eliminate the free evolution of equation ( p3[ ) , before taking the mean value and making use of 
Wick's theorem and the RPA to find the equation of motion 

^^ d{alan) ^ 2Y,clc,{xm\V\an){n^ - n„)e-^"--*. (47) 

The solution of equation (^7|) differs slightly from that of the anomalous and linear terms we considered earlier. We 
previously made the assumption that the approximate solutions were zero at time i = — c», but this is not the case 
here. The diagonal term (to = n) of (al^an) is the average number of particles in mode n. For a system at finite 
temperature in which the incoherent region begins with some population, this will have a non-zero initial value, and 
appears in the solution of equation (B^) as a constant of integration. We find 

(ai^an) = -27ri^c*Co-(x"i|F|crn)(n,„ - n„)6{hujxman) + n„iS„in, (48) 

where in the spatial representation of the FTGPE the term involving the Kronecker delta function represents the 
mean field of the incoherent region acting on the coherent region wave function. We have neglected the principal part 
in the solution as it is dominated by the energy conserving processes. 



In the basis set representation of the FTGPE the corresponding term ( 28f ) is therefore 



dc 
in^ = ... + 2j2nmCp{am\V\mP) (49) 

771/3 

- Awi ^ C/3(Q;TO|y|n/3)^c*c<x(xm|y|(Tn)(n„ ~ nn)6{huj^^aa) 

mnj3 X"' 

The processes it describes are analogous to the scattering term that was introduced in the model of condensate growth 
in references [ pO| , pl| . It represents an incoherent particle colliding with a coherent particle, with the coherent particle 
moving between levels within the region C . In the presence of a condensate this process is recognized as Landau 
damping, and in references pQ,pli it was shown to have an important effect in the onset of condensate growth. 
However in simulations closer to equilibrium when a condensate is already present, the off-diagonal part of this term 
is unlikely to be large as the forward and backward rates will be similar. 

D. The grow^th term 



Finally, we consider the term involving {iffifj) on the fourth line of the FTGPE (29d) which wc identify as the 
growth term. This will generally be the most important term involving the field f) in the FTGPE, as it will be 
responsible for the majority of particle exchange between the coherent and incoherent regions. While the linear 
terms and the anomalous terms can be important in certain circumstances near T = 0, in most situations at finite 
temperature they are small in comparison with the growth term because of the size difference between the coherent 
and incoherent regions. 

To find an explicit expression for this term we begin with the equation of motion for aXcimCim equation ( |C4| ). 
Eliminating the free evolution, taking the mean value and making use of Wick's theorem and the RPA leaves us with 

d{aXaman) v-^ 

ifi ^^ ==2 2^ cx{mn\V\qx)e''^"^^''-^^{nq{l + n,„ + n„) - n„,nn}. (50) 

X 

As this can represent energy conserving processes, the approximate solution is 

{a\a,nan) = 2'Ki'^c^{mn\V\qx)5{TiuJmnqx){nr,inn -n,(l + n,„ +n„)}. (51) 



Substituting equation ( |51| ) into the basis set version of the FTGPE we find the term corresponding to equation ( |28g 
becomes 
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ih —^ = ... + 2TTi ^ c^{aq\V\mn){mn\V\qx)S{h,uJmnqx) 

qmnx 

x{n„n„ -nq(l + n„ + n„)}. (52) 

We can understand the physical meaning of this result by comparing it with the equivalent quantum Boltzmann rate 



for the growth of a coherent region level with population N^ (see equation A14). Making the approximation that 

(l + Na) ~ Na, we find 



" ex ^(1 + Na){l + nq)nmnn - Nang{l + nm){l +nn), 

qmn 



dt 



^ 7Vc<[n„n„ -nq(l + n„j + n„)], (53) 



and we can see that the right-hand sides of (|5^) and ( p3| ) are similar. 

Equation (JS^) is very similar to the growth part of the description of condensate formation of reference [pO| . 
However, rather than describing the rate of change of an occupation number of a quasiparticle level, it describes the 
growth in amplitude of a basis component making up the coherent region wave function i/i. To calculate this term for 
inclusion in the FTGPE requires both a reasonably good basis and corresponding energies for the coherent region, 
along with a method of calculating or approximating the matrix elements that appear in equation (|5^) . While this 
is not difficult in principle, in practice they need to be calculated at each time step in the evolution of the FTGPE. 
This issue will be treated in future numerical work. 

VI. THE PROJECTED GROSS-PITAEVSKII EQUATION 

The finite temperature Gross-Pitaevskii equation contains all the necessary elements for a complete description of a 
condensed Bose gas, given that the occupation number conditions are satisfied. However, it is somewhat complicated 
to implement numerically. Some insight can be gained by neglecting all the terms coupling the coherent region to 
the effective bath 77, and considering the equation of motion for the coherent region region alone — the first line of 
equation (p9|). We call this the projected Gross-Pitaevskii equation (PGPE) 

'^^^ = H,,i,{^) + UoV {|V(x)P^(x)} , (54) 

and we have studied solutions of this in reference ||l^. Although it is not immediately evident, the PGPE does 
conserve normalization and energy, and this can be easily understood by considering the effective Hamiltonian as we 
discuss further in appendix pi 

The projected GPE describes a microcanonical system, whereas in the full FTGPE the region C will fluctuate in 
energy and number of particles. However, if the region C contains sufficiently many modes, then fluctuations in energy 
and particle number in the grand canonical ensemble would be small. Hence we would expect an equilibrium state of 
the projected GPE to be similar to that of the finite temperature GPE coupled to a bath f)(x) with the appropriate 
chemical potential and temperature. 

The projected GPE by itself cannot capture the entire physics of the Bose field at finite temperature. Indeed, in 
any system in which there is a significant thermalized coherent region that may be modelled by this equation, there 
will be a much larger incoherent region whose effects will be important. Nonetheless, we showed in reference p^ ] 
that the GPE without any additional terms can describe evolution of general configurations of the coherent region C 
towards an equilibrium that can be parameterized by a temperature. 

The detailed non-equilibrium dynamics of the system will depend on the exchange of energy and particles between 
C and the bath. However, we expect that for modelling many experiments at finite temperature the projected GPE 
alone may be sufficient. It is well known from kinetic theory that particles mainly interact with others of similar 
energy. If the presence of the bath is important, then it could be well approximated by a constant temperature and 
chemical potential, alongthe lines of the formalism outlined in reference pO| and used in the model of condensate 
growth in references [pO|,pl|. This issue will be considered further elsewhere. 
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VII. CONCLUSIONS AND PROSPECTS 

In this paper we have derived an approximate formahsni for calculating the dynamics of a thermal Bose gas in the 
highly-occupied limit. We have derived a non-perturbative finite temperature Gross-Pitaevskii equation that describes 
the evolution of the wave function of the coherent region C, and identified and discussed the physical meaning of each 
of the terms that arise. In particular we have indicated how the anomalous term introduces the T-matrix in coherent 
region collisions, and that this can describe elastic particle loss in condensate collisions at large relative momenta. 



The terms analogous to the scattering and growth processes of reference 30 1 have also been discussed. 

While the formalism discussed here leaves out the possibility of quantum correlation effects in BECs, it should 
nevertheless provide a basis for the description of many experimentally accessible features that can not be described 
by the ordinary Gross-Pitaevskii equation. An important aim of our development has been to produce a numerically 
tractable formalism, and the first step in a practical implementation has been made in reference |l7| where it was 
shown that the projected GPE will evolve a generalized random initial distribution to an equilibrium described by a 
temperature. The next step is to include some of the additional terms of the FTGPE described in section M. 

We can also expect the FTGPE to describe the phase transition region, as long as the condition Nk ^ 1 is satisfied 
for the coherent region modes. The physics of phase transitions is generally classical in nature, being dominated by 
large fluctuations at long wavelengths. This of course is what the GPE describes, and in fact it has been used as a 
model of phase transitions in other areas of condensed matter physics. The GPE has the same energy functional as 
that used in the classical renormalization group theory of the superfluid phase transition, and it is therefore reasonable 
to expect the same approximations to be valid for the case of BEG. Our formalism, which couples the GPE to a kinetic 
treatment of the thermal particles, provides a non-perturbative dynamical finite temperature theory of BEG. This 
can then be used to study the region of the phase transition where perturbative approaches fail, and indeed where no 
other techniques are available. 
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APPENDIX A: DERIVATION OF THE QUANTUM BOLTZMANN EQUATION 

In this appendix we give a derivation of the quantum Boltzmann equation, which gives an accurate description 
of the time evolution of a Bose gas well above the transition temperature. In this regime the mean time between 
particle collisions is long compared to the duration of a collision, so the eigenstates of Hq provide a good basis and 
the interaction part of the Hamiltonian Hj can be treated as a perturbation. This is the method by which incoherent 
region terms are treated in this paper. 

The operators dp have no mean value above the transition temperature, and so we want an equation of motion for 
the mean number of particles in mode p, (hp) — {dj,dp). From equation (Um wc find 

^ = - J E {pq\V\mnyaldldradne'<^'^''+'^^-^--^-^' + h.c, (Al) 

qmn 

where h.c. is the hermitian conjugate. In order to find a closed expression for the evolution of fip, we begin by finding 



an equation of motion for the quantity dj^d^d^dn which appears on the right-hand side of equation (Al). However, this 
new equation for four operators contains terms involving six operators, and the equations of motion for six operators 
involve eight operators and so on. We proceed by truncating this series at the first iteration, approximately solving 
the equation for dj^d'ljdmdn, and substituting the result back into equation (Al). 



One way to derive an equation of motion for djjd'^dmdn would be to commute it with the Hamiltonian. An equivalent 
method (which will be useful for other purposes later) is simply to use the chain rule. We can formally write the 
solution as 



ifln = / '^^''jp (°P^ 



dldldmdn = / dt'-^ (dldldmdn) 



15 



dt' 



dal 



— n'n n 






-—a^^aman + 0-}p-7—araan + % 



dt 



dt' 






dt' 



~i ~i ~ dan 
+ alalam^ 



(A2) 



We now substitute equation dlQ) for each of the ddk/dt. This is a straight-forward but tedious process, which we 
illustrate only on the last term of equation (A2). We have 



dt I a^ajjam—- 

jkl 



ij.' ~i ~i ~ ~t~ ~ i(uJi+uj^—uju—Ljj)t' 

dt a^pttljamajakaie ^ *^ ^ " , 



dt' 



--}_^S,n{^J\V\kl) 
jkl 

x{dldld^jdrndkdi + dld\akdi5mj)t 



i{uli+ujj—uJk—u>i)t' 



(A3) 



where we have arranged the creation and annihilation operators in normal order. 

We now introduce our first approximation. The free evolution of the operators has already been removed and so 
assuming the interaction is a perturbation, over the period of the integral most of the time dependence is contained 
in the exponential. Invoking the Markov approximation (which assumes that correlations between the operators are 
not important on the time scale of interest), we can therefore take the operators outside of the integral and replace 
their time dependence by the time at the upper limit. The integral then leaves us with a delta function when the 
frequencies of the modes add to zero, and a principal part when they do not — a standard result normally expressed 
as 



-J^ = p(l] -iTrS(x). 
x + ie \xj ' ' 

We assume that the principal part is negligible, which leaves us with 



(A4) 



dt' 



' ' "\0'\am-Tp- ) = T- E 5in{ij\V\kl)5{u^ + LUj - CJfe - Wi) 



jkl 



x{dld\d]dradkdi + a],ajafea;(5„j). 



(A5) 



Combining all the terms that arise out of this procedure, and making use of the symmetries of the indices, equation ( Al) 
becomes 



ijkl 



UJl) 



dla'jdldidmdnSkq + dldjdldidmdnSkp + dldjdmdnSkqSip 
-dldldldmdkdidjn - dldldldndkdiSjm ~ a)pd\dkdi5jn5v 



(A6) 



The next step is to take the expectation value of equation (A6), as we are concerned with the time evolution of 
the average occupation of any individual level. In fact, it is most useful to take the ensemble average, which can be 
written as 



{A) = TripA), 



(A7) 



where p is the density matrix of the system, and Tr denotes the trace operation. We are left to calculate quantities 
such as {dldjdmdn)- To do so, we assume that our system is sufficiently near thermal equilibrium that we can use 
Wick's theorem ||32|1 . This states that for any system with a Hamiltonian that is a quadratic form in creation and 
annihilation operators, the ensemble average of any product of these is simply the contraction of all possible pairings 
of the operators. For example we have 
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{a\a}-aman) = {a\a}j) {amCin) + {a\an) {a^^am) + {a\am) {a]ar, 



(A8) 



This relation is exact at thermal equilibrium, and should be a very good approximation nearby. 

We now make use of a further approximation, known as the random phase approximation (RPA) , which states that 
if the density matrix for the system is diagonal then 



{a\a]) = {h,hj) = 0, {a\hj) = UiSij. 



(A9) 



This will be an excellent approximation when there is no condensate present, as the interaction Hamiltonian Hj is 
only a small perturbation to the ideal gas Hamiltonian Hq. Thus we have 






On substituting these relations into equation (A6) the final result is 



-^ = ^Yl \{P<l\y\mn)fS{ujp+ujq 

qmn 



UJn) 



Stt 



X [{np + l){ng + l)nmn„ - npnq{nm + !)(«« + 1)] 
'Y{pq\V\pn) {mn\V\mq)npnra{nn - nq)S{uJn ~ t^q) 



(AlO) 
(All) 



(Al2a) 
(A12b) 



qmn 



(Al2c) 



The first p art of this expression (Al2a) is the standard quantum Boltzmann equation; however, the lines of (A12b) 
and (Al2c) do not appear in most definitions of the QBE. We would like to note the following about these terms: 



1. The scattering processes described by the matrix elements in th ese ter ms involves a third particle, and hence 
these collision terms are of higher order than those described by (Al2a). 



2. If we calculate the matrix elements using the contact potential approximation in the homogeneous limit, then 
these terms become 



8t:U, 



(Al2b) -^ 2o2 X! ^^'^1 ^ k„)^npn„(n„ - nq)6{u!n - ^Uq), 



(A12c 



8ttU§ 



y^ <5(kp - \<infnqnm{nn - np)6{ljJn - UJp). 



(Al3a) 
(A13b) 



The delta functions in momentum are equivalent to Kronecker delta functions in the quantum labels for the 
system, and hence these terms vanish. 

3. For an ergodic system where the occupation of a level depends only on its energy, these terms are again identically 



4. The delta functions in frequency depend on only two of the partic le indi ces, r ather than four as for equa- 
tion (Al2a). This means there will be far fewer matches for ( Al2b| ) and ( Al2c ), and therefore these can be 
considered surface terms that become small in the thermodynamic limit. 

We are therefore justified in neglecting these terms, and are left with the usual quantmn Boltzmann equation 



-if- = ^ X! l(wl^l™'^)P'^('^p + <^9 -'^™ -^") 



qmn 



x<{np + l){nq + l)n„in„ - UpUq^rim + l)("-n + 1) 



(A14) 
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Validity 

We summarize the validity conditions for the QBE as: 

1. The Markov approximation must be vahd such that correlations induced by collisions are unimportant. 

2. The system should close enough to equilibrium such that the factorization of Wick's theorem is valid. 

3. There must be a good basis such that the RPA is valid. In our derivation wc have assumed that Hj should be 
a perturbation to the system for this to hold. However, if the average effect of Hj can be absorbed into Hq to 
form an effective Hamiltonian with a good basis, then the QBE derivation may still be valid. It must be noted, 
however, that the Markov approximation may not be valid in this regime as the mean collision time will be 
much reduced. 

APPENDIX B: THE GPE KINETIC EQUATION 

It is interesting to consider the kinetic equation that would result if we assume that the GPE is a good description 
of the system of interest. We expand the time dependent wave function as 

V'(x,t)=^5fc0fc(x)e-''^''*, (Bl) 

k 

where the {ck} are slowly varying. Substituting equation (Bl) into the time-dependent GPE (m and performing the 
operation J d'^x0*(x) on both sides results in the basis set representation of the GPE 

in^ = ^(M|y|mn)a;£„5„e^(--+-'"-" -"")*. (B2) 

qmn 



We note that equation (B2) is identical in form to the basis set equation of motion for the Bose field (16) but for 



the replacement Uk <-^ Ck- In fact, equation (B2) could be derived directly from equation (|l^) by expanding the 
annihilation operators as ak = Ck + 5k and then neglecting the quantum fluctuation terms 5k ■ 



Wc can now carry out the same procedure on (B2) as was applied to equation ( pq ) in the derivation of the QBE. 
The only difference is that we are now manipulating c-numbers rather than operators, and so any terms arising from 
commutators in the previous treatment will not appear. This means that the terms of the form d^a^aa will disappear 
from equation ( |A6| ), leaving only terms involving six c-numbers. Writing Up for c*Cp, the resulting GPE kinetic 
equation is 

dup 47r/7o v^ I , I \i2<r/ , N 

-j^ = -^— 2_^ \{pq\mn)\ d{ep + e^ - e„i - e„) 

qmn 

X <^ {Up + nqjrijnrin - npnq{n,n + Un) \ (B3) 



which is exactly the same form as the QBE (A14) except that the spontaneous collision terms are excluded. This 
equation was first considered by Svistunov in a study of the formation of a condensate in a weakly-interacting Bose 
gasi. 

Some of the approximations made in the derivation of the GPE kinetic equation may not hold in the presence of 
a condensate. In particular, the assumption of correlations being unimportant on the scale of the collision time is 
unlikely to be valid. The GPE kinetic equation does, however, give us an understanding of the collision processes 
that are described by the full GPE. 



From equation (B3) we can see that the GPE contains stimulated collision processes only. To understand this. 



consider the collision p -I- q ^ m -I- n. This process will only be represented by the GPE if one of the levels (to, n) is 



already occupied. This is in contrast to the QBE, for which the term in the curly brackets of equation (A14) can be 
written 

(1 + Up -I- nq)nmnn - npnq{l + n„i + n„) >. (B4) 

Thus, due to the neglect of the quantum nature of the modes, the GPE can only accurately describe the evolution 
and interaction of modes which satisfy Up ^ 1, such that (1 + np + nq) k, [up -f Uq). 
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dt 



APPENDIX C: INCOHERENT REGION EQUATIONS 

In this appendix we write out the Heisenberg equations of motion for all the operator combinations that appear in 
the basis set version of the FTGPE (p8|). The single operator term is 

(Cla) 
(Clb) 
(Clc) 

(Cld) 

(Cle) 

(Clf) 

(Clg) 



+ {pip\V\'ii)ip) 

+ ^{pq\V\il;^)al 

+ 2^(j}^\V\mtl))a,r, 

m 

+ y^(p-0|l^|™n)a„ia„ 
+ 2y^ij)q\V\mil))a\a^, 

qm 

+ '^{pq\V\mn)a\am 

qmn 



The equation of motion for aj, is simply the hermitian conjugate of (CI). The other equations of motion can be found 
either by calculating the commutator with the Hamiltonian, or using the chain rule. We have 



in = n(UJjn + LUnjamari 



dt 



+ {rmjj\V\ip'ip)an + {nip\V\->p'ip)am 
+ {mn\V\'ipip) 

+ 2J {kn\V\ip')p)alam + {km\V\'ip^)alar, 
k 

+ 22J [{nij}\V\k'4!)aka„i + {mil)\V\kilj)akan] 
k 

+ 51 [("^V'|V'|fcj}a„afcaj + {nil)\V\kj)a,nakaj] 

kj 

+ 2> {mn\V\kip)ak 

k 

+ 2 2, {km\V\jip)alajan + {kn\V\jilj)al.aja„ 
kj 

+ y ^{■mn\V\kj)akaj 

kj 

+ ^ [{mq\V\kj)alanakaj + (nq|F|fcj}aJa,„afeaj] 

qkj 



(C2a) 

(C2b) 
(C2c) 

(C2d) 
(C2e) 
(C2f) 
(C2g) 
(C2h) 
(C2i) 
(C2j) 



in = n{UJn - OJmjalnan 



dt 



+ {m(j\V Itp ^p)al,^ — {^pip\V\m^p)dn 

+ ^ {kn\V\ijJ'il;)aldl^ - {tpil;\V\km)dkar, 

k 

+ 22J {'ml^\V\ktlj)al^dk - {k2lj\V\rmlj)dlan 

k 

+ ^ \{ni:\V\kj)a'l^dkaj - (fcj|y|m-0)a|.a]a„ 
kj 



(C3a) 
(C3b) 
(C3c) 

(C3d) 

(C3e) 
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+ 



2^ ^{kn\V\ji})alnalaj - {ji)\V\km)a^.aka. 
y^ \{qn\V\ki)al^a\akaj - {kj\V\qm)a},ja\aqan 



(C3f) 
(C3g) 



and 



d{alaynan) 
in 2_^ = in 



dt 



^^ d{d,ndn) d{a\) ,^ ^ 



h{u!m + tOn - UJq)a)qdmdn 

{mi!\V\'4!%l))d\dn + {n^\V\i!tj;)d\dm 
{mn\V\'ip^)d\ - {il;^\V\qi^)dmdn 

k 

- {i>il^\V\kq)dkdyndn] 
2^ [(nVl^|fc'0)aJafeam + {m'ip\V\kij) d\dkdr, 

k 

- {kTp\V\q'il))d\dmdn 
y^ [(mV'|F|fcj)aJa„afeaj + (nV'|F|fcj)atamafcaj 

-(fcj|y|gV)a|a]a„ia„ 

2 > {mn\V\k%l])a)dk 
k 

2 ^ (fcm|V^|j-0)aJaJ,aja„ + (fcn|T/^|j-0)aJaJ,aja„ 

-(jV'|V^|fcg)a^%a„ia„ 
y^(mn|T/|fcj)a^afcaj- 

kj 

y^ [(mr|F|fcj)aJata„afcaj + (nr|y|A:j)aJa];a,„afeaj 
-(fcj|F|rq)aj^^a]a,.a,„a„ 



(C4a) 
(C4b) 
(C4c) 

(C4d) 

(C4e) 
(C4f) 

(C4g) 
(C4h) 

(C4i) 
(C4j) 

(C4k) 

(C41) 
(C4m) 

(C4n) 

(C4o) 



APPENDIX D: INCOHERENT REGION KINETIC EQUATION 

In this section we present the kinetic equation of motion that can be derived for the incoherent region, using the 
same methods as in the derivation of the QBE in appendix ^ and the techniques outhned in section M. We find 



r = -T-^\^<l\y\0-P)? K^^pqap)\CaCl3\^ ijlp + Uq + I) 



dt 



qaf3 



Stt 



-^ X! l(P'^l^l'^/3)l^'^(^^pam/3)|CaC/3|^(n™ -Tip) 



af3m 



47r 



y^ \{pa\V\'mn)fS{hLJparnn) 
amn 

x|cQp{n,„n„ ~np{l + n,„ +n„)} 



(Dla) 
(Dlb) 

(Die) 
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+ -^ TJ \{pq\V\ma)\^5{hujpqrna) 



aqm 

2 



X |cq| {{rip + rig + l)nm - ripUq} (Did) 



47r 



, y^ \{pq\V\mn)\'^S{hLUpq,nn) 

qmn 

X {{1 + np){l + ng)nmnn ~ ripng{l + nm){l + Tin)} ■ (Die) 

We can recognize each of these terms from cm" previous discussions: 



1. The term (Dla) is from the anomalous term, and is only non-zero when we consider the collision of multiple 
condensates. 



2. The line (Dlb) describes the scattering processes. 



3. The two terms (Die) and (Did) are due to the forward and backward growth terms of the coherent region. 



4. The line (Die) is simply the QBE for the incoherent region. 



It may seem surprising that there is no contribution from the linear terms discussed in the previous section. This is 
because the rates for each of the forward and backward processes contain only stimulated terms and so they cancel. 
The kinetic equation for the incoherent region is thus the usual QBE but with additional couplings to the coherent 
region whose physical meanings can be understood from the corresponding terms in section |V|. 

APPENDIX E: CONSERVATION OF NORMALIZATION IN THE PGPE 

The PGPE conserves normalization and energy because the effective projected Hamiltonian is hermitian. The 
nonlinear term of the GPE can be considered to describe interactions between two particles, and as such there can 
be collisions in which two coherent atoms collide and one is ejected from the coherent region. However, the projector 
excludes these terms from the equation of motion which we now demonstrate. 

Consider the equation of motion in a basis set. By substituting rpi^} — X^feec '''='?^fe (■''■) ™^'-' equation ( |54|) and 
performing the operation J d'^X(/)p(x) on both sides we find 

ih—^ = hujpCp + UqP ^ c*gC,nCn(jjq\mn). (El) 

qmn^C 

If the state p € C then all four labels are from the coherent region and there is no transfer of population outside the 
region. For p ^ C the matrix element {pq\rnn) is not zero necessarily, and therefore it seems collisions between states 
from within the coherent region can transfer population outside of C. However, we should not be considering the 
equations of motion for amplitudes p ^ C in the first place, as they not in the definition of the wave function ijj{x). 
The projection operation is therefore performed implicitly by the basis set representation. 

Numerically solving the GPE using a basis set method requires a triple summation over indices, which is a very 
time-consuming operation. This suggests that we should instead use the spatial representation of equation (p4), where 
the nonlinear term is local. However, any spatial grid that is fine enough to provide a good representation of all the 



states within C will also be able to represent modes outside the region C. From equation (El) we can see that this 
will cause population to be transferred outside of C, and so in this case we need to consider the projection operation 
explicitly. 

Another method of approaching this issue is to assume that all modes in the problem can be represented by the 
GPE, but artificially choose part of the system to be the coherent region. Thus the field operator can be written as 

'l'(x)«^(x) + r;(x), (E2) 

where both fields are classical. Substituting this into the interaction part of the Hamiltonian (|5|) and using the contact 
potential approximation, we have 

Hi/Uo = H4 + H3 + H2 + Hi+Ho, (E3) 

Hi = ^m^i^, (E4) 
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TABLE I. Classical FTGPE divided into terms representing physical processes involving n coherent states 
coherent ift(^j=Px... ih(^\=Qx 



No. of 

;oheren1 

states 



"4 \i^ 

3 +2|V|^?? + r;>2 +1^1^ 

1 +\v\^V +2V'|»7P + ^Z'*'?" 

+0 +\riW 



i/g = xjj*'ilj*ipri + 'iJj*r)*iP^, (E5) 

-^2 = ni^* ip* W + '2^'* T]* ''Pv + TiV* V* '^'4' y (E6) 

i/i = i]j*rj*riri + rj^rf^r], (E7) 

^0 = ^r;*r?*,7,7, (E8) 

where we have dropped all space and time labels for clarity, and have divided the Hamiltonian into five terms depending 
on the number of coherent region fields they contain. Considering the Hamiltonian in this form we can easily interpret 
each of the terms. Each '^p* creates a particle in the coherent region, and each ij] removes a particle from the coherent 
region. The rf and 77 perform the same operation outside the region C . This allows us to identify which processes 
each of the terms in the Hamiltonian contribute to the equations of motion for both ^ and 77. 
We can now derive equations of motion for ^ and 77 by using functional differentiation. We find 

As an example, let us consider the contribution to these equa tion s for all interactions involving three coherent and 



one incoherent state described by H^. We find from (|E5|) and (E9) 



in^^r{2\^\^v + v*i^^), (eio) 

^n^^Q{\iJ\'^J). (Ell) 

The results of carrying out this operation for all particle processes are summarized in table ffl. The equations of 
motion for the system will together conserve both energy and normalization if all terms in any row of the table are 
included, as this correctly accounts for all forward and backward processes of the same order. We have confirmed this 
numerically by carrying out coupled simulations of ^p and rj and including only some of these terms. 

We can now see that if we want an equation describing interactions involving four coherent states, but neglecting 
all processes involving incoherent particles, then this is simply the PGPE. 
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